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(February 7, 2008) 

Abstract 

Deep underground muon events recorded by the Soudan 2 detector, located at 
a depth of 2100 meters of water equivalent, have been used to infer the nuclear 
composition of cosmic rays in the "knee" region of the cosmic ray energy 
spectrum. The observed muon multiplicity distribution favors a composition 
model with a substantial proton content in the energy region 8 x 10 5 — 1.3 x 10 
GeV/nucleus. 
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I. INTRODUCTION 



The composition of the cosmic rays in the "knee" (~ 10 4 TeV/nucleus) region of the 
cosmic ray all-particle spectrum has consequences for astrophysical models of particle ac- 
celeration and propagation. For example, a model of accretion onto a black hole in the 
center of an Active Galactic Nucleus predicts an excess of protons around the knee region 
in the cosmic ray flux |JJ, while a model of shock acceleration during a supernova blast into 
a stellar wind environment predicts an excess of heavy nuclei in the cosmic ray flux in the 
same energy region M. 

Unfortunately, the composition at energies above ~ 1000 TeV/nucleus is difficult to 
measure directly due to the steeply falling spectrum of cosmic ray primaries. The flux of 
particles with energies greater than 1000 TeV/nucleus is only ~ 60/m 2 /sr/year. Therefore 
measurements of cosmic rays in this energy regime require detectors of either very large 
acceptance or long exposure, neither of which is currently feasible for measurements near 
the top of the atmosphere or above. Instead, studies of the composition of cosmic ray 
primaries in this high energy regime can be carried out indirectly through measurements 
of aspects of the atmospheric cascade generated by the interaction of a cosmic ray primary 
with the earth's atmosphere. 

Deep underground experiments, such as Soudan 2, indirectly study the composition of 
cosmic rays by comparing observations of multiple muon event rates to expectations derived 
through Monte Carlo simulations using various trial composition models as input. A multiple 
muon event is one in which two or more nearly parallel, time-coincident muon tracks are 
observed in the detector. These muons are decay products of mesons which are generated in 
the hadronic core of the atmospheric cascade. At high energies, massive primaries generate 
more high energy muons per event than proton primaries of the same total energy. This 
is because the initial parent meson (predominantly pion) particle multiplicities are larger 
from the more massive primary, and the atmosphere is more favorable to pion decay versus 
interaction early in the cascade development. In addition, the point of first interaction of 
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the heavy primary is likely to be higher in the atmosphere than for the proton primary due 
to the larger cross-section of the heavy nucleus-air interactions. 

At the Soudan site, the measurement of the high energy muon component underground 
is coupled with sampling of the electromagnetic component of the air shower at the surface 
using a small proportional tube array |J , and a Cherenkov light detector array [Q . A correct 
interpretation of the cosmic ray composition should yield consistent results using all three 
experimental techniques in all possible combinations. 

In this paper, we report on an analysis in which the observed multiple muon event rates 
recorded in the Soudan 2 detector are compared with the simulated rates obtained using 
three distinct trial composition models. The next section of this paper contains a brief 



description of the relevant aspects of the detector. In Sec. Ill, we report on the analysis of 



the observed multiple muon event rates. Sec. [iy| contains a discussion of the Monte Carlo 
simulation. Sec. [V] has a discussion of the test composition models used in this analysis 
and a comparison between the data multiple muon rates and the rates observed using these 
models. Finally, Sec. [VT] summarizes the results. 



II. THE DETECTOR 

The Soudan 2 detector || is a high resolution tracking calorimeter located in the Soudan 
Underground Mine State Park in northern Minnesota, at a depth of 710 meters below the 
surface of the earth. This depth corresponds to a muon threshold energy of ~ 0.7 TeV for 
a muon transmission probability of at least 50%. The modular design of the detector has 
allowed the continuous acquisition of data from the beginning of its construction in mid- 1988 
to the present, during which time Soudan 2 has recorded more than 33 million muon events. 
The detector reached its full size of 224 1 mxl mx2.7 m modules in November, 1993, for a 
final operating size of 8 mxl4 m in horizontal surface area x5.4 m in height. 

The active region of each detector module consists of 7560 15 mm diameter plastic drift 
tubes layered between 241 1.6 mm thick corrugated steel sheets. Ionization deposited in the 
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tubes drifts under the influence of an electric field toward the tube ends, where it is collected 
on vertical anode wires spaced 15 mm apart, and horizontal cathode strips spaced 10 mm 
apart. The pulses on the wires and strips are read-out and digitized every 200 ns, which, for 
a drift rate of 0.6 cm//is, corresponds to a timing resolution along the drift axis of 1.2 mm. 
(The measured resolution along the drift-axis of ~ 6 mm is larger than the timing resolution 
due to variations in drift velocities.) Pulse matching of the anode and cathode signals yields 
a three-dimensional space coordinate for each drift tube crossing along the path of a charged 
particle. 

Events with high muon multiplicity require the ability to separate tracks bunched tightly 
together, as well as to distinguish tracks obscured by showers in multiple muon events in 
which a large bremsstrahlung has occured. The high resolution of Soudan 2 is particularly 
suited to this type of study since a typical track will leave hundreds of reconstructed pulses 
along its path. Fig. p] shows three views of a 14-muon event as recorded in the detector. 
The resolution of a space point along a muon track detected in Soudan 2 is ~ 1 cm. The 
angular resolution of a typical muon track is < 1°. 

The primary trigger operating over the time span of this data sample was called the 
"edge" trigger. It required a muon to have pulsed a minimum of 7 anode wires or 8 cathode 
strips out of any contiguous block of 16 channels of either type, separated by at least one 
trigger clock pulse of length 600 ns, and occuring within the given time window of 72/xs. 
The effect of the trigger requirement on a muon track of length 50 cm, the minimum track 
length considered in this analysis, is shown in Fig. |2]. The boundary shown in the figure 
represents the ideal; in actuality, muon bremsstrahlung or pair production initiated showers 
assist triggers allowing some muons outside of this region to satisfy the trigger requirements 
as well. This, combined with individual channel inefficiencies, creates some fuzziness at the 
trigger acceptance boundaries. To avoid the trigger holes, the data events in this analysis 
were subjected to software imposed zenith and azimuthal angle cuts as discussed in the next 
section. 

In addition to the main detector, there is a 14mx31mxl0m veto shield consisting of 
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proportional tubes which nearly covers the entire An steradian surrounding the main detec- 
tor. For multiple muon analysis, the shield is useful as a tool for selecting multiple muon 
event candidates in the main detector. It can also be used to study muon tracks which pass 
outside of the main detector volume ||, however this has not been included in the analysis 
presented here. 

III. DATA ACQUISITION AND ANALYSIS 

The data set reported in this paper consists of a sample of the total number of events 
recorded in Soudan 2 over the time span June 1991 through October 1991, and includes 
7.2xl0 5 single muon and 22000 multiple muon events after all cuts are applied. The detector 
size at the time of this data sample was 8 mxll m in surface area x5.4 m high. A larger 
data sample is not necessary for this analysis, since its accuracy is dominated by systematic 
uncertainties in the Monte Carlo simulation and not by statistics. 

Data in Soudan 2 are accumulated in modular "runs" each lasting a little over an hour 
and containing typically 1500-2500 events. To be included as part of the data sample, each 
run needed to pass a series of "run quality" checks in order to exclude runs with localized 
hardware failures such as high voltage trips, excessive noise, or anomalously high or low 
average trigger rates. After application of run cuts, the total live-time considered in this 
analysis was 1348.8 hours. 

Several cuts were applied to the data sample to ensure the highest quality data. Each 
muon track was required to have a minimum length of 50 cm. Multiple muon events were 
required to satisfy a "parallel" cut such that for each muon track there was an angular 
separation of less than 5° from at least one other muon in the group. This restriction was 
to eliminate locally produced tracks. To avoid the trigger holes, we have made cuts on 
the azimuthal and zenith angle regions of acceptance. The muon events are confined to the 
azimuthal regions 10° < <p < 80° in each quadrant and zenith angles > 15°. The muon events 
are further confined to zenith angles < 60° as this allows us to apply the flat atmosphere 
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model which is prevalent in atmospheric cascade simulations. The zenith and azimuthal 
angles of an event with more than one muon track were determined as the average of all 
muon tracks satisfying the length and parallel cuts. 

By design, the Soudan 2 track reconstruction software finds and reconstructs everything 
from through-going muons to "contained" tracks consisting of as few as 5 pulses. This 
is because Soudan 2's primary purpose is to search for the short tracks left by proton 
decay. The software is very efficient at finding and reconstructing single muon events which 
typically contain several hundred pulses. For this reason, determining the single muon event 
rate required only purifying the total of all tracks found by the reconstruction software to 
extract the through-going muon tracks. The regular track reconstruction code was modified 
to tag muon tracks as those which satisfied at least one of the following criteria on both ends 
of the track: 

• The track extrapolated to a time-coincident veto shield hit. 

• The last reconstructed hit on the end of the track was within 50 cm of the detector 
edge. 

• The end of the track projected through a detector crack, which is defined as the small 
open region between each pair of modules. 

These requirements were enough to make the software very efficient at producing a very 
pure single muon sample. The single muon event reconstruction was tested against 1900 
randomly selected events from two separate data runs, of which 594 events were found to be 
single muon events passing the length and angular cuts considered here. The results of this 
test are shown in Table |. The software was found to be 99.2 ± 0.4% efficient at identifying 
single muon events which passed the stated angular and length cuts. The background 
of misreconstructed tracks contributing to the single muon sample was determined to be 
2.0 ± 0.6%. These corrections have been applied to the observed number of single muon 
events in Table piT[ 
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The software reconstruction of multiple muon events in the main detector is complicated 
due to a hardware design feature which electrically adds signals from separate regions of 
the detector together during readout. The purpose of this "multiplexing" is to decrease 
the cost of the electronics required to read out the 4 x 10 6 drift tubes, and it occurs just 
before the pulses are digitized. The multiplexing is configured so that a match between 
a given anode pulse and cathode pulse has a unique position in the detector. The pulses 
are demultiplexed at the software stage. This multiplexing has little effect on the software 
reconstruction of contained events or single muon events because the relatively small number 
of pulses in these events leads to a simple interpretation of the data, but as the multiplicity 
of the event increases, the number of pulses and the complexity of the demultiplexing of the 
event increases as well. For this reason, the software reconstruction of all candidate multiple 
muon events was supplemented by scanning performed by a physicist using an interactive 
graphics program. 

Multiple muon candidates were selected based on a set of generous but simple criteria. 
These criteria were that a candidate event contain both of the following: 

• At least one "good" track satisfying the criteria of the single muon tracks described 
above, except that in this case the angular restrictions on this one track were loosened 
to 12.5° < 9 < 62.5° and 7.5° < < 82.5° in each quadrant. 

• At least one 2-dimensional anode-time or cathode-time reconstructed track, which, 
when paired with the opposite 2-d track from the "good" track, was parallel to the 
"good" track within 5°. The sum of the projected lengths of all parallel 2-d tracks in 
either the anode-time or cathode-time view had to be at least 50 cm. 

The net effect of these cuts was to select ~ 6% of the total number of muon events as multiple 
muon candidates, while the final sample post-scanning consisted of only 3% multiple muons. 
Therefore these selection cuts are considered to be very liberal. 

The multiple muon candidates were examined not just for multiplicity but to check their 
reconstruction results which were displayed directly over the event. The majority of low 
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multiplicity events (N^ < 4) had been reconstructed correctly by the software, while the 
highest multiplicity events (> 8) generally required some manual corrections which were 
applied through the interactive graphics program to produce the correct fits. 

Events with multiplicities greater than 12 were discarded in this analysis due to a flaw 
in the data acquisition hardware which existed at the time of the data sample. This flaw 
put a maximum limit on the amount of readout time taken by any one event from the 
time of the trigger to the end of the event readout CAMAC sequence. The amount of 
time spent in this sequence rises with the size of the event, which is in turn dependent 
on the event multiplicity. Since the readout time distributions corresponding to a given 
multiplicity have long tails, there is a gradual lessening in the efficiency to readout the 
entire event with increasing multiplicity. Through a simulation in which low multiplicity 
data events (for which the readout time distributions are known) were used to generate the 
readout time distributions of high multiplicity events, the efficiency for reading out high 
multiplicity events was determined 0. These efficiencies were calculated based on a worst 
case scenario, and as such cannot be used to correct the data. In this simulation, 99% of 
12 muon events were read out successfully by the data acquisition software. This is the 
maximum multiplicity considered in this analysis. 

The efficiency of the multiple muon candidate selection criteria is tabulated in Table 
0. In this case, the same two runs used to determine the single muon efficiency were used 
to determine the number of multiple muon events which satisfy the angular, length, and 
parallel cuts described above. Of the 2970 events in these two runs scanned for multiple 
muon events, 30 events were determined to be multiple muon events which satisfied these 
cuts, and the multiple muon selection software found all 30 of these events. Five additional 
runs shown in the table were scanned by a team of undergraduate scanners. Seventy-nine out 
of the 8830 events scanned were found to be multiple muon events which satisfied the cuts, 
of which the multiple muon selection software found 79. The efficiency of the multiple muon 
selection software has therefore been determined in excess of 97.5% at the 90% confidence 



level. This uncertainty has been applied to the errors in Table |TJ for all multiplicities, even 



though the selection software will certainly be more efficient at finding high multiplicity 
events than low multiplicity events, so that this can be considered to be a very conservative 
estimate of the total uncertainty at high multiplicities. 

The final corrected muon event rates observed in Soudan 2 are shown in Table [TTI| and 
Fig. §. 

IV. MONTE CARLO SIMULATION 

The increasing statistical accuracy of indirect measurement data correlated to phenom- 
ena from primary energies in the knee region requires the use of sophisticated Monte Carlo 
simulations to constrain the systematic errors in these types of measurements. In this ap- 
plication, we have used the most fully-developed Monte Carlo simulation of the atmospheric 
cascade and muon rock propagation currently available. The Monte Carlo simulation used 
in this analysis consists of three stages: 

• A 3-dimensional simulation of the atmospheric cascade. 

• A 3-dimensional simulation of the propagation of the muons through the rock over- 
burden. 

• A simulation of the detector. 

The cascade simulation uses the HEMAS cascade code || for controlling the overall 
structure of the cascade development. This code injects a nucleus of a requested energy, 
mass, azimuth and zenith angle into the atmosphere. It then tracks this particle and any 
secondaries produced along the path of the cascade development until they either interact, 
decay, drop below some user defined energy threshold, or reach the atmospheric sampling 
height. In our case, this last quantity is at the surface of the earth above Soudan 2, which 
corresponds to 492 meters above sea level. 

The HEMAS cascade code is built so that the hadronic interaction code is fully modu- 
lar. In this analysis, we have used the recent program SIBYLL |§ for generating hadronic 
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interactions. SIBYLL is based on the dual parton model with minijet production superim- 
posed. It was designed with the motivation of using a theoretical model for extrapolation 
of "low" energy accelerator and fixed-target data to the energies necessary for the study 
of cosmic-ray interactions. SIBYLL agrees reasonably well with available accelerator data 
with one notable exception: the (p?) distributions associated with particle production in pp 
interactions are underestimated at large Feynman x (x^p > 0.15) [f|. The (pr) distributions 
of parent mesons are significant in studies like this because (pt) is the dominant contribut- 
ing factor to the lateral distribution of underground muons. At the Soudan 2 depth, the 
mean separation of muon pairs is comparable to the size of the detector. This means that 
the Monte Carlo simulation of muon lateral spread is significant for a detector of our size 
and depth since we need to determine correctly the rate of "observed" events of a given 
multiplicity from the total number of events at our depth. We have tested the effect that 
the SIBYLL underestimation of (pr) has on our analysis and discuss this in Sec. |V|. 

The nucleus-air interaction simulation is also a modular component of the HEMAS cas- 
cade code. We have used the NUCLIB ||TU[ nuclear interaction code for this stage. In 
addition, the HEMAS cascade simulation package has been modified from its original form 
to include the effect of the local geomagnetic field, which has a strength of 0.59 Gauss and a 
magnetic inclination and declination of 75.1°N and 0.85°E respectively at the Soudan 2 site. 
Even though the earth's magnetic field is weak, it plays a noticeable role in determining 
the transverse displacement of the underground muons because the distance traveled by the 
cosmic ray muons through the atmosphere is very long. Because the magnetic field near the 
Soudan site is nearly vertical, particles at large zenith angles are more greatly affected than 
those near vertical. We have found that the addition of the magnetic field to the cascade 
simulation has a negligible effect for events with zenith angle 9 = 15°. However, at 9 = 60°, 
the earth's magnetic field increases the mean transverse displacement of muons from the 
event core at the Soudan 2 depth from 13 m with no magnetic field to 18 m with magnetic 
field. 
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To determine the absolute single and multiple muon rates in the detector, it is essential 
to understand the composition and thickness of the rock overburden, as well as to have a 
means to simulate the passage of the muons through the rock. A digitized map from the 
US Geological Survey was used to determine the depth of rock around the Soudan site. 
The density of the rock at the Soudan site is also variable. Soudan is located in the "Iron 
Range" of northern Minnesota and the rock overburden is of non-uniform composition with 
pockets of iron ore interspersed among an overburden consisting mostly of Greenstone JIT 



To determine the density of the overburden, a fit to "world survey" muon data [E| was 
performed using a large body of muon data spanning the periods June, 1991 to March, 
1996. The effective rock density was determined in 337 angular bins covering the angular 
region considered in this analysis 0. 

The GEANT Detector Description Simulation |13[ package developed by CERN was used 
to propagate muons underground. Muon energy loss mechanisms are fairly well understood 



TM , and in fact have been experimentally verified up to energies of ~ 1 TeV for muon energy 



loss in iron [15,16"! . For Soudan 2 data this is the critical region of energy, since it defines the 
shape of the rise of the transmission probability curve. (The transmission probability curve 
is the probability of a muon to successfully reach the Soudan 2 level versus muon surface 
energy.) We have compared GEANT to the available muon experimental energy loss data 
with good results |L7| . 

Finally, the detector was simulated using a simplified model which compares well against 
a more realistic simulation of the response of the tracking calorimeter modules and their 
readout electronics. All cuts applied to the data were also applied to the Monte Carlo 
simulated events. 



V. ANALYSIS 
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A. Composition Models 



The energy and mass of the primary cosmic rays cannot be determined on an event-by- 
event basis through indirect experimental methods such as underground muon rates because 
of large fluctuations in the atmospheric development of a cosmic ray cascade. Instead, the 
experimental approach is to assume a model of primary composition at the top of the 
atmosphere as a function of mass and energy, and to use the Monte Carlo simulation to 
predict from the assumed composition the muon rate underground. For simplicity, primary 
composition models generally divide primary cosmic rays into five mass groups centered 
around the principal nuclei H, He, CNO, Ne-S, and Fe. To narrow the field of possible test 
composition models, we have defined three composition models which satisfy the following 
criteria: 

• The model is theoretically motivated by an astrophysical model. 

• The model fits the available satellite and balloon direct measurement data in the low 
energy region (< 1000 TeV/nucleus) in which this data is available. 

• The model normalizes to the air-shower determined all-particle spectrum in the knee 
energy region. 

In recent years the amount of direct measurement cosmic ray composition data, both 
satellite and balloon, has grown considerably. Direct measurements currently extend up to 
almost 1000 TeV/nucleus ji~8,19||. Monte Carlo simulations show that Soudan 2 should be 



sensitive to muon events generated by primary cosmic rays in the energy region 5 - 50000 
TeV/nucleus. Therefore, there is significant overlap between direct measurements and the 
energy region under study in our analysis. 

Use of the available direct measurement data performs two important functions in our 
analysis. First, it tightly constrains the possible test composition models, e.g. a composi- 
tion model of pure protons over the entire energy range under study would be in obvious 
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contradiction to the mixed composition observed in the lower energy regime by direct mea- 
surements. Secondly, proper normalization to the direct measurement data allows for a test 
of the atmospheric cascade simulation at low energy. To illustrate, Fig. |] shows the fraction 
of events which are produced by primaries of energy less than 100 TeV/nucleus as a function 
of multiplicity. It is clear from the figure that greater than 85% of single muon events are 
generated by primaries in an energy region in which the composition is well known. There- 
fore, the absolute single muon rate in Soudan 2 can be used as an important test of the 
Monte Carlo simulation. 

As has been pointed out elsewhere ||20|| , the availability of new high energy direct mea- 



surement composition data has made some popular composition models somewhat obsolete. 
In particular, we note that the versions of a Light test composition model used by NUSEX 
||2T|| and MACRO [^] no longer give good agreement to the available direct measurement 
data over the entire relevant energy range. To formulate test composition models, we follow 
the lead of Silberberg, et al. |23|] and Stanev, Biermann, and Gaisser |Q and assume a 



two-component model at low-energies, such that the differential spectrum is described by 

dN 

— = K X E~^ + K 2 E-^ (1) 
dhi 

for each of the five mass groups. A two component model follows naturally from the as- 
sumption that the low energy cosmic rays come primarily from two sources: supernova blasts 
into a homogeneous interstellar medium(e.g. ||25|| ) and supernova blasts into a stellar wind 
environment |p26[. The latter is theorized to play a significant role for the heavy elements 
and to produce a flatter spectrum than that of the former |2^j23| . To fit Eq. (JJ) to the direct 



measurement data, we have compiled data from a large number of direct measurement ex- 
periments. These measurements are shown in Figs. |5| and []. We found that good agreement 
between Eq. ([[]) and all five mass groups could be obtained by using 7x=2.75 and 7 2 =2.50, 
and by allowing the coefficients Ki and K 2 to have the values tabulated in Table [IV]. Eq. 
(0), along with the values in Table [TV], is used to describe the low-energy component in all 
models considered in our analysis. The agreement of this low-energy component with the 
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direct measurement data is shown in Figs. [| and || 

We then define three test composition models: "New Source_P" , "Heavy" , and "New 
Source_Fe". The New Source_P model, motivated by the model of Fichtel and Linsley 
f2?|| , assumes that a new source consisting entirely of protons predominates at energies 
around the knee. This type of proton-rich source is compatible with the ideas of the AGN 
particle acceleration model of Szabo and Protheroe Q. In the New Source_P model, we have 



assumed that the low-energy components have an exponential cutoff (as in [24J]), such that 
the differential spectrum for these components becomes 

^ = {K-yE-^ + K 2 E-^)e~^tt. (2) 
dE 

The exponential cutoff energy, E cni , is determined for the proton component by the direct 
measurement data. The exponential cutoff for the heavier elements is allowed to extend 



out to higher energies, as suggested by reference j26| for a stellar wind component, and as 
required by the helium data. This cutoff for the heavier elements is also taken to be rigidity 
dependent. The low-energy cutoff for all five mass groups are given in Table [V] and are 
shown in Fig. ||. 

The New Source_P high-energy component, shown in Fig. |7|, has the functional form 

^ = K e~ A ° E ~ B ° E~ 7 °, (3) 



with parameters for this component given in Table VI. We have chosen to normalize the 



high-energy component such that the summed mass components of the model equal the all- 



particle spectrum as determined by Akeno |[28|| . This choice of normalization and its effect 



on the underground muon rates will be discussed in the next section. 

The second test model, the Heavy model, is motivated by the theoretical model of Bier- 
mann, et al. P,|24|j, in which the stellar wind component extends out to energies up to 
~Zx 10 8 GeV/nucleus. In the Heavy test model used here, each of the heavy mass groups 
(He,CNO,Ne-S,Fe) has the low energy form given in Eq. (|IJ) up to a sharp break at energy 
Efc nee , above which the mass groups follow the single power- law 
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— oc E ». (4) 



The Efc nee and 73 parameters for this model are given in Table |V11| and the agreement of this 



model with direct measurements is shown in Fig. |6|. The normalization to the all-particle 
spectrum closely matches that of New Source_P and is shown in Fig. |7|. 

The third test model, New Source_Fe, is exactly the same as New Source_P except that 
the high energy protons have been replaced by iron to test the opposite extreme. 

The fractional composition for all three models is shown in Fig. || We note that even 
though a pure proton high energy component is assumed for the New Source_P model, this 
model is still fairly heavy in what will be the critical energy region for this study: 10 6 — 10 7 
GeV/nucleus. In particular, we note that the New Source_P model is heavier than the 
MACRO Light |22| model in this energy region. 

B. Results 

The results of this analysis are shown in Fig. [5| This figure shows the ratio of simulated 
multiple muon event rates to observed event rates for each of the three test composition 
models. Fig. |lTj shows the primary energy range to which each multiplicity is sensitive. We 
note that 90% of the high multiplicity (> 6) events come from a fairly well defined region 
of the energy spectrum which is 8 x 10 5 — 1.3 x 10 7 GeV/nucleus. This is the energy region 
just below and around the knee of the energy spectrum. From Fig. ^|, therefore, it is clear 
that the high multiplicity events offer a statistically significant test of composition of cosmic 
rays in this part of the all-particle spectrum. 

The agreement of the measured and simulated event rates of each test composition model 
is tabulated in Table |V11J|. x 2 /d.o.f. in each case is calculated by comparing the absolute 



event rates of the simulated and measured data for the seven multiplicities 6 through 12 
which span the primary energy region of interest around the knee of the all-particle spectrum. 
The data event rates used in this calculation include statistical and systematic errors, while 
the simulated event rates include statistical errors only. Of the three trial composition models 
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considered here, the New Source_P simulated event rates clearly give the best agreement with 
the data in the knee energy region of the all-particle spectrum. 

The observed absolute single muon rate agrees with the Monte Carlo rates obtained 
with all three composition models to within 5%. As already mentioned, ~ 85% of single 
muon events come from a region of the spectrum in which the composition is known from 
direct measurements; therefore the agreement of the single muon observed and Monte Carlo 
rates shows that the Monte Carlo is successful at simulating muon events at low energy. Of 
course, as we extend the measurement to higher energy, the systematic uncertainties become 
more complex. Here we discuss two of the more important contributors to the systematic 
uncertainty of the Monte Carlo simulation at energies around the knee of the all-particle 
spectrum. 



1. All-particle spectrum normalization 
As already mentioned, we have chosen to normalize the test composition models to 



the Akeno determined all-particle spectrum [28| in the knee energy region. The choice of 



normalization in this energy range has some arbitrariness associated with it since the all- 
particle spectrum is determined indirectly through air-shower measurements. This choice 
obviously has an effect on the high multiplicity underground muon rates. The Akeno group 
f28fl established a range in which the all-particle spectrum might fall based on their own 
measurements using various techniques and the measurements of other experiments. Their 
parametrized form of their measured all-particle spectrum, used in the analysis presented 
here, lies along the low end of this range. If we were to increase the normalization of the 
test composition models towards the upper end of their range, the effect would be to push 
the results further away from the heavy composition models. 

We also compare the normalization of the test model all-particle spectra to the recently 
reported result of the Tibet AS7 [Q air shower collaboration. The Tibet air-shower array 
operates at the high altitude of 4300 meters above sea level, corresponding to an atmospheric 
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depth of 606 g/cm 2 . This is an ideal altitude for studying the all-particle spectrum in the 
energy region just before and around the knee, since air showers generated by primaries in 
the energy range ~ 10 5 — 10 8 GeV reach their maximum shower development close to the 
altitude of the experiment, and hence their fluctuations at the point of sampling 
minimum. Also, the average shower sizes generated by primaries of different atomic mass 
become nearly identical at sampling heights near shower maximum. This is very important 
for eliminating the systematic dependence of the measurement of primary energy on the 
assumed primary composition. 

The all-particle spectra of the test composition models considered in this analysis lie 
about 22% lower than the Tibet measurements at energies of 10 5 ' 75 GeV, and about 9% 
higher at energies of 10 7 GeV. A renormalization of the all-particle spectrum of the Heavy 
model to the Tibet spectrum yields underground high-multiplicity event rates which are 
~ 10% higher than the present rates, again pushing the Heavy model further away from the 
data. 

2. Hadronic interaction model 

As already noted, one of the uncertainties in the hadronic interaction model is that 
SIBYLL tends to underestimate the increase in (pp) with xp (xp > 0.15) for particle pro- 
duction in pp interactions. On the other hand, an older version of a Monte Carlo hadronic 
interaction code, HEMAS ||, produces a significantly greater mean transverse momentum 
in the forward fragmentation region than that of SIBYLL, and is in fact likely to overes- 
timate the (pp) at large x^ 0. We have run a test using the Heavy model as input and 
SIBYLL as the main driving code, but scaling the transverse momentum in each SIBYLL 
hadronic interaction according to the HEMAS (pp) distributions. In this way, we were able 
to extract and test the effect that just this one aspect of the hadronic interaction code has 
on the Monte Carlo simulated Heavy model rates. The result is that approximately 10% 
fewer Monte Carlo events with multiplicities > 6 are seen in the detector given the larger 
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(jPt) distributions of the HEMAS hadronic interaction code. This lowers the Heavy compo- 
sition model rates in the direction of the data, resulting in an improved x 2 /d.o.f. agreement 
between the data and the Heavy model. The x 2 /d.o.f. for the (pr) scaled Heavy model rates 
(calculated, as before, by comparing the absolute event rates of the simulated and measured 
data for the seven multiplicities 6 through 12) is 55/7 (C.L. = 1.8xlCT 7 %). 



VI. CONCLUSION 

We have shown that the multiple muon rates observed in the Soudan 2 detector are 
sensitive to the nuclear composition of cosmic rays in the energy region 8 x 10 5 — 1.3 x 10 7 
GeV/nucleus, which is the energy region just before and around the knee in the cosmic ray 
all-particle spectrum. The composition model favored in this work includes an enhanced 
component of protons in this energy region. This component is compatible with the AGN 
model of particle acceleration of Szabo and Protheroe |1[], which shows a contribution of 
protons due to AGN sources localized in the knee energy region of the spectrum. 

The light composition model favored in this analysis is compatible with previous results 
from Soudan 1 and Soudan 2 [Q, in which a small surface array was operated in con- 
junction with the underground detector. The result is also in agreement with the recent 



multiple muon analysis performed by MACRO |45| , as is shown in Fig. [Ty. Results from 
Fly's Eye |30| favor a mixed to heavy composition at energies near 3 x 10 8 GeV/nucleus. We 
note that the Fly's Eye result does not directly contradict our result, since our measurement 
applies to lower primary energies. 
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FIG. 1. A 14 muon event in the Soudan 2 detector. The figure shows three different view of 
the same event. All axes are labeled in centimeters. 
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FIG. 2. The shaded regions correspond to the "idealized" trigger holes of the detector as seen 
by muons with a track length of 50 cm. The coordinate system is defined such that ((f) = 0°,9 = 90°) 
points North, ((f) = 90°, 6 = 90°) points West, and 9 = 0° points straight up from the detector floor 
plane. 
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FIG. 3. The absolute event rate as a function of multiplicity as seen in Soudan 2. The rates 



also tabulated in Table TT|. 
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FIG. 4. The fraction of events at each multiplicity generated by primaries with energy less 
than 100 TeV/nucleus. 



27 



Proton 



Helium 





1.0 2.0 3.0 4.0 5.0 6.0 

CNO 



.0 2.0 3.0 4.0 5.0 6.0 

Ne-S 




1.0 2.0 3.0 4.0 5.0 6.0 

Fe Group 



1.0 2.0 3.0 4.0 5.0 6.0 

Log 1 0(Energy(GeV/nucleon)) 



F 



1.0 2.0 3.0 4.0 5.0 6.0 



Log 1 0(Energy(GeV/nucleon)) 

FIG. 5. Differential fluxes of the five mass groups used in the low energy component of the New 
Source_P and New Source_Fe composition models. Muon events observed in Soudan 2 are products 
of primaries with energies above ~ 10 3 GeV/nucleon; however, it is useful to show differential flux 
measurements at energies below this to display trends in the data. Below ~10 GeV/nucleon, the 
observed flux of cosmic rays is greatly affected by solar modulation, and the amplitude of the data 
will depend upon the intensity of solar modulation at the time the measurement was taken. The 
fits for each of the five mass groups, shown as solid curves, have been "modulated" at low energy 
in this figure by a factor of the form (l+aE _,s ) _1 to fit the data, however this does not affect the 
low energy parametrizations given in the text for the energy region > 10 3 GeV/nucleon which is 
of interest to us. The direct measurement data have been compiled from the following sources: o 
[31]; • p^] ; □ |33| for proton and helium and [34] for heavier elements; filled □ [35|; A |36| for 



proton and helium and [37] for heavier elements; filled A [19|; V | IS ] ; filled y 38]; O p9[; filled O 



|40fl; * fl4lH; and X |]42|. The "Fe Group" consists of the elements Z=26-28 for o and □; Z=26 for 
filled A; and Z=26-30 for A. 
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FIG. 6. Differential fluxes of the five mass groups used in the Heavy composition model. The 
references for the direct measurement data can be found in the caption for Figure [|. 
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FIG. 7. The New Source_P and Heavy composition models as compared to the all-particle 
spectrum. The all-particle spectra of the two composition models have been normalized to the 
Akeno parametrized all-particle spectrum |28[ (hashed region) in the knee region. The rest of the 
all-particle data are from the compilation by Stanev |43(| . The five mass groups shown are for P 
(bold), He (dash), CNO (dot), Ne-S (dot-dash), and Fe (solid). The New Source_Fe model is the 
same as New Source_P with the high-energy proton component replaced by iron. 
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FIG. 8. The fractional composition of the three test composition models as a function of 
energy. 
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FIG. 9. Ratio of simulated to data absolute event rates as a function of event multiplicity. 
The ratio errors are calculated using the systematic and statistical errors of the Data rates, and 
statistical errors only of the Monte Carlo rates. See the text for a discussion of Monte Carlo 
systematic errors. 
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FIG. 11. Average atomic mass versus primary energy for each of the three test composition 



models, and as determined by MACRO [45]. 
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TABLES 



TABLE I. Efficiency and background for reconstructed single muon events. 



Run 


Date Total Scanned Muon Events 
Events Passing Cuts 


Efficiency 




Background 


29052 
31111 


7/21/91 
10/14/91 


1000 273 
900 321 


§i = 0.993 ± 0.005 
§±f = 0.991 ±0.005 


4 
275 

8 
326 


= 0.015 ±0.007 
= 0.025 ± 0.009 


Total 




1900 594 


§| = 0.992 ± 0.004 


12 
601 


= 0.020 ± 0.006 


TABLE II. Efficiency of the software at finding multiple muon events satisfyin 
ified in the text. 


g the cuts spec- 


Run 


Date 


Total Events 


Multiple Muon Events 


Efficiency 


29052 
31111 
28132 
33009 
35018 
35520 
35848 


7/21/91 
10/14/91 
6/13/91 
1/01/92 
3/31/92 
4/26/92 
5/12/92 


2070 
900 
1858 
1868 
1795 
1428 
1881 


23 
7 
12 
13 
19 
16 
19 




1. 
1. 
1. 
1. 
1. 
1. 
1. 


Total 




11800 


109 




109 i 
109 ~~ ±- 
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TABLE III. Number of events and event rate observed at each multiplicity. 



Multi- 


Events 


Events ± stat ± syst 


Rate (hr x )± stat ± syst 


plicity 


(uncorrected) 


(corrected) 


(corrected) 


1 


724792 


716024 ± 846 ± 5249 


530.87 ± 0.63 ± 3.89 


2 


17237 


17237 ± 131 ± 417 


12.780 ± 0.097 1[J- 309 


3 


2813 


2813 ± 53 tf 


2.086 ± 0.039 to 050 


4 


905 


905 ± 30 ±f 


0.671 ± 0.022 t° - 016 


5 


385 


385 ± 20 l^- 4 


0.285 ± 0.015 t° m7 


6 


172 


172 ± 13 to 1 


0.128 ± 0.010 t° 003 


7 


102 


102 ± 10 tl 5 


0.0756 ± 0.0074 t° 0019 


8 


58 


58 ± 7.6 +l A 


0.0430 ± 0.0056 t 0010 


9 


42 


42 ± 6.5 


0.0311 ± 0.0048 +0- 0007 


10 


37 


37 ± 6.1 t°o 9 ° 


0.0274 ± 0.0045 t° 0007 


11 


20 


20 ± 4.5 ±°' 48 


0.0148 ± 0.0033 +° 000A 


12 


13 


13 ± 3.6 t° 31 


0.0096 ± 0.0027 +° 0003 



TABLE IV. Low-energy parameters of the New Source_P, New Source_Fe, and Heavy compo- 
sition models. The units of K±,2 are [m _2 s _1 sr~ 1 (GeV/nucleus) 71 ' 2 ^ 1 ]. 



Mass Group 


Zeff 


Aeff 


Ki 


7i 


K 2 


72 


H 


1 


1 


20830. 


2.75 


0. 


2.5 


He 


2 


4 


7750. 


2.75 


840. 


2.5 


CNO 


7 


14 


3545. 


2.75 


550. 


2.5 


Ne-S 


12 


24 


2655. 


2.75 


445. 


2.5 


Fe 


26 


56 


2120. 


2.75 


335. 


2.5 
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TABLE V. The low-energy component energy cutoff used in the New Source_P and New 



Source_Fe composition models. E cut is 


given in (GeV/nucleus). 


Mass Group 




E cu t 


n 




1.5 x 10 5 


He 




1.0 x 10 6 


CNO 




3.5 x 10 6 


Ne-S 




6.0 x 10 6 


Fe 




1.3 x 10 7 


TABLE VI. 


High-energy parameters of the New Source_P and New Source_Fe composition 


model. The units of K G are [m 2 s 1 sx 


_1 (GeV/nucleus) 7 " l ] and A G is given in (GeV/nucleus) B °. 


K 


A 


B 7o 


2.75E9 


327. 


0.322 3.3 


TABLE VII. 


High-energy parameters of the Heavy composition model. Efc raee is given in 


GeV/nucleus. 






Mass Group 




Efcnee 73 


He 




1.2 x 10 6 3.08 


CNO 




4.2 x 10 6 3.08 


Ne-S 




7.2 x 10 6 3.08 


Fe 




15.6 x 10 6 3.08 
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TABLE VIII. Calculated values of x 2 , X 2 P er degree of freedom (d.o.f.), and confidence level 
for each of the test composition models. See text for further discussion. 

Composition Model x 2 d.o.f. X 2 / d -°- f - C.L. (%) 

New Source_P 8.7 7 1.2 28 

Heavy 103 7 15 2.1 xl(T 17 

New Source_Fe 236 7 34 3.0xl0" 45 
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